//============================================================================
// Name        : Coevolution 1
// Author      : Alex Seeholzer
// Description :
// Benedikt Obermayers Original Code, catalytic System of 4 cyclically catalyzing wild types
// including a functional region
//============================================================================

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>

#define TYPE unsigned int
#define L 32
#define ONE 1
#define RRMASK 0x000000FF
#define MAXINT 0xFFFFFFFF
#define LRR 8
#define DELTA 1.
#define N 10000
#define TMAX 1000

#define NM 4 /* population number */
const TYPE masters[] = {0x00000000, 0xFFFFFF00, 0xAAAAAA00, 0x55555500,
			0x11111100};

const double alpha=5.;
const double gam=NM*10.;

inline TYPE mutate (TYPE x, double mu);
inline int hd (TYPE x, TYPE y); // hamming dist
inline void outbinary(TYPE num);
void reproduce (TYPE *pop, double mu, unsigned int cdis[NM][(ONE << LRR)]);
TYPE compute_consensus (TYPE *pop);

// computes for a pop the distribution of clones (stored in cdis)
void compute_cdis (TYPE *pop, unsigned int *cdis, TYPE master);

// Random number generator
gsl_rng *r;
gsl_interp_accel *acc;
const gsl_rng_type *T;

int main(int argc, char **argv) {

	// Set up output file
	FILE * outFile;
	outFile = fopen("out.dat","w");

	int i,j,k,t,s;
	TYPE pop[N]; // array of population
	double mu; // mutation rate
	TYPE smax[NM]; // maximal populated sequences for each population
	unsigned int cdis[NM][(ONE << LRR)]; // for each master, distribution of functional region

	// random number generators
	gsl_rng_env_setup();
	acc = gsl_interp_accel_alloc();
	T = gsl_rng_gfsr4;
	r = gsl_rng_alloc(T);

	// read mu, s = seed
	if (argc > 2) {
		mu=atof(argv[1]);
		s=atoi(argv[2]);
		fprintf(stderr, "mu=% .3e, s=%i\n", mu, s);
	} else {
		fprintf(stderr, "no mutation rate given!\nUsing mu=% .3e, s=%i\n",mu=.0005,s=23);
		//exit(1);
	}

	gsl_rng_set (r, s); // initialize random num generator with seed s

	// population initialization
	for (j = 0; j < NM; j++)
		for (i = j * N / NM; i < (j + 1) * N / NM; i++)
			pop[i] = (masters[j] & ~RRMASK)
					| gsl_rng_uniform_int(r, ONE << LRR);

	// compute frequency of master sequences for initial reproduce call
	for (j = 0; j < NM; j++)
		compute_cdis(pop, cdis[j], masters[j]);

	// main loop - reproduce, calc distribution and output
	for (t = 0; t < TMAX; t++) {
		reproduce(pop, mu, cdis);
		for (j = 0; j < NM; j++)
			compute_cdis(pop, cdis[j], masters[j]);

		for (k = 0; k < (ONE << LRR); k++) {
			fprintf(outFile,"%i %i ", t, k);
			for (j = 0; j < NM; j++)
				fprintf(outFile,"%i ", cdis[j][k]);
			fprintf(outFile,"\n");
		}
		fprintf(outFile,"\n");
		if (t%100 == 0)
			printf("time: %i\n",t);
	}


	for (j = 0; j < NM; j++) {
		smax[j] = 0;
		for (k = 1; k < (ONE << LRR); k++)
			if (cdis[j][k] > cdis[j][smax[j]])
				smax[j] = k;

		fprintf(stderr, "smax[%i]= %3i (%02X), neighbors: ", j, smax[j],
				smax[j]);
		for (k = 0; k < LRR; k++)
			fprintf(stderr, "%i ", smax[j] ^ (ONE << k));
		fprintf(stderr, "\n");
	}

	fclose(outFile);
	return 0;
}

void compute_cdis (TYPE *pop, unsigned int *cdis, TYPE master) {
  for (int k=0; k<(ONE << LRR); k++)
    cdis[k]=0;

  for (int i=0; i<N; i++)
    if ((pop[i] & ~RRMASK)==(master & ~RRMASK))
      cdis[(unsigned int)(pop[i] & RRMASK)]++;

  return;
}

void reproduce(TYPE* pop, double mu, unsigned int cdis[NM][(ONE << LRR)]) {
	int i, j, index;
	TYPE popnew[N];
	double ptot;
	double sel[N];
	double mut[N];
	double select[2*N];
	double delta[NM]; // wieviele anteilig in pop j den genotyp von typ i haben (in RR) == Xi

	for (i = 0; i < N; i++) {
		for (j = 0; j < NM; j++)
			delta[j] = cdis[j][(pop[i]) & RRMASK] / (double) N;

		mut[i] = 1.;
		sel[i] = 1/DELTA; /* 1./DT */

		for (j = 0; j < NM; j++)
			if ((pop[i] & ~RRMASK) == (masters[j] & ~RRMASK))
				mut[i] = alpha + gam * delta[(j + 1) % NM];
	}

	/* sum of probabilities for normalization */
	ptot = 0;
	for (i = 0; i < N; i++)
		ptot += sel[i] + mut[i];

	/* evaluate probabilities to be selected
	 (even index) and to be selected+mutated
	 (odd indices) */
	select[0] = sel[0] / ptot;
	select[1] = mut[0] / ptot + select[0];
	for (i = 1; i < N; i++) {
		select[2* i ] = sel[i] / ptot + select[2* i - 1];
		select[2* i + 1] = mut[i] / ptot + select[2* i ];
	}
	/* get N random numbers in [0,1) */
	for (i = 0; i < N; i++) {

		/* look up entry in select array */
		index = gsl_interp_bsearch(select, gsl_rng_uniform(r), 0, 2*N - 2) + 1;

		/* if index is odd, mutate this sequence with
		 probability mu, otherwise copy  */
		popnew[i] = (index % 2 == 1 ? mutate(pop[index / 2], mu) : pop[index
				/ 2]);
	}

	for (i = 0; i < N; i++)
		pop[i] = popnew[i];

	return;
}

inline TYPE mutate (TYPE x, double mu) {
  int i;
  TYPE out=x;
  for (i=0; i<L; i++)
    if (gsl_rng_uniform(r) < mu)
      out^=(ONE << i);
  return out;
}

inline int hd (TYPE x, TYPE y) {
  int dist=0;
  TYPE val=x ^ y;
  while (val) {
    ++dist;
    val &= val-1;
  }
  return dist;
}

inline void outbinary(TYPE num) {
	TYPE tmp = num;
	for(int i=0;i<L;i++) {
		int powr = int(pow(2,(L-i-1)));
		int res = tmp/powr;
		fprintf(stderr,"%i",res);
		tmp-=res*powr;
	}
	fprintf(stderr,"\n");
}
